#include "Cell.h"
//#include "ParallelData.h"
#include "Point.h"
#include "Facet.h"

using namespace dolfin;

//-----------------------------------------------------------------------------
double Facet::normal(uint i) const
{
  const uint D = _mesh->topology().dim();
  _mesh->init(D - 1);
  _mesh->init(D - 1, D);
//  dolfin_assert(_mesh->ordered());

  // Get cell to which face belong (first cell when there is more than one)
  const Cell cell(*_mesh, this->entities(D)[0]);

  // Get local index of facet with respect to the cell
  const uint local_facet = cell.index(*this);

  return cell.normal(local_facet, i);
}
//-----------------------------------------------------------------------------
Point Facet::normal() const
{
  const uint D = _mesh->topology().dim();
  _mesh->init(D - 1);
  _mesh->init(D - 1, D);
//  dolfin_assert(_mesh->ordered());

  // Get cell to which face belong (first cell when there is more than one)
  const Cell cell(*_mesh, this->entities(D)[0]);

  // Get local index of facet with respect to the cell
  const uint local_facet = cell.index(*this);

  return cell.normal(local_facet);
}
//-----------------------------------------------------------------------------
bool Facet::exterior() const
{
//  if (_mesh->parallel_data().exterior_facet().size() > 0)
//    return _mesh->parallel_data().exterior_facet()[this->index()];
//  else
    return num_entities(dim() + 1) == 1;
}
//-----------------------------------------------------------------------------
std::pair<const Cell, const Cell>
Facet::adjacent_cells(const MeshFunction<uint>* facet_orientation) const
{
//  dolfin_assert(num_entities(dim() + 1) == 2);

  // Get cell indices
  const uint D = dim() + 1;
  const uint c0 = entities(D)[0];
  const uint c1 = entities(D)[1];

  // Normal ordering
  if (!facet_orientation || (*facet_orientation)[*this] == c0)
    return std::make_pair(Cell(mesh(), c0), Cell(mesh(), c1));

  // Sanity check
  if ((*facet_orientation)[*this] != c1)
  {
    printf("Facet.cpp: extract adjacent cells of facet. Illegal facet orientation specified, cell %d is not a neighbor of facet %d",  (*facet_orientation)[*this], index());
  }

  // Opposite ordering
  return std::make_pair(Cell(mesh(), c1), Cell(mesh(), c0));
}
//-----------------------------------------------------------------------------
